Efficient and accurate three dimensional Poisson solver for surface problems 
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We present a method that gives highly accurate electrostatic potentials for systems where we have 
periodic boundary conditions in two spatial directions but free boundary conditions in the third 
direction. These boundary conditions are needed for all kind of surface problems. Our method 
has an 0(N log N) computational cost, where N is the number of grid points, with a very small 
prefactor. This Poisson solver is primarily intended for real space methods where the charge 
density and the potential are given on a uniform grid. 
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I. INTRODUCTION 

Electrostatic potentials play a fundamental role in 
nearly any field of physics and chemistry. Having ef- 
ficient algorithms to find the electrostatic potential V 
arising from a charge distribution p or, in other words, 
to solve the Poisson's equation 



V 2 V = -Aitp , 



(1) 



is therefore essential. The large variety of situations in 
which this equation can be found lead us to face this 
problem with different choices of the boundary conditions 
(BC). The long-range behavior of the inverse Laplacian 
operator make this problem to be strongly dependent on 
the BC of the system. 

The most immediate approach to the Poisson equation 
can be achieved for periodic BC, where a traditional re- 
ciprocal space treatment is both rapid and simple, since 
the Laplacian matrix is diagonal in a plane wave rep- 
resentation. If the density p is originally given in real 
space, a first Fast Fourier Transformation (FFT) is used 
to transform the real space data in reciprocal space. The 
Poisson equation is then solved in reciprocal space and 
finally the result is transformed back into real space by a 
second FFT. Because of the FFT's, the overall computa- 
tional scaling is (0(N log N)) with respect to the number 
of grid points N. 

The situation is different if one considers the same 
problem for free BC. In this case the solution of Pois- 
son's equation can formally be obtained from a three- 
dimensional integral: 



V(v) 



J dr'Gflr- 



■r'l)p(r'), 



(2) 



where G(r) = 1/r is the Green function of the Lapla- 
cian operator in the unconstrained R 3 space. The long 
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range nature of the kernel operator G does not allow us to 
mimic free BC with a very large periodic volume. Con- 
sequently, the description of non-periodic systems with 
a periodic formalism always introduces long-range inter- 
actions between super-cells that falsify the results. Due 
to the simplicity of the plane wave methods, various at- 
tempts have been made to generalize the reciprocal space 
approach to free BC (@; O; [Til ). All of them use a 
FFT at some point, and have thus a 0(N log N) scaling. 
These methods have some restrictions and cannot be used 
blindly. For example, the method by Fusti-Molnar and 
Pulay is efficient only for spherical geometries, and the 
method by Martina and Tuckerman requires artificially 
large simulations boxes that are expensive numerically. 
Nonetheless, the usefulness of reciprocal space methods 
has been demonstrated for a variety of applications, and 
plane-wave based codes are widely used in the chemical 
physics community. 

Another choice of the BC that is of great interest is for 
systems that are periodically replicated in two dimen- 
sions but with finite extent in the third, namely surface 
systems. The surface-specific experimental techniques 
developed in recent years produced important results {]]), 
that can benefit from theoretical prediction and analysis. 
The development of efficient techniques for systems with 
such boundary conditions thus became of great impor- 
tance. Explicit Poisson solvers have been developed in 
this framework [j; la), with a reciprocal space based 
treatment. Essentially, these Poisson solvers are built 
following a suitable generalization for surfaces BC of the 
same methods that were developed for isolated systems. 
As for the free BC case, screening functions are present to 
subtract the artificial interaction between the super-cells 
in the non-periodic direction. Therefore, they exhibit the 
same kind of intrinsic limitations, as for example a good 
accuracy only in the bulk of the computational region, 
with the consequent need for artificially large simulation 
boxes which may increase the computational overhead. 

Electrostatic potentials can either be calculated by 
solving the differential Poisson equation or by solving 
the equivalent integral equation cq.([2]). The methods 
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that solve the differential equation are iterative and they 
require various tuning. A good representative of these 
methods is the multigrid approach (J7|). Several different 
steps such as smoothing, restriction and prolongation are 
needed in this approach. Each of these steps has to be 
tuned to optimize speed and accuracy. Approaches based 
on the integral equation are in contrast straightforward 
and do not require such tuning. 

In this paper we will describe a new Poisson solver com- 
patible with the boundary conditions of surfaces. Con- 
trary to Poisson solvers based on reciprocal space treat- 
ment, the fundamental operations of this Poisson solver 
are based on a mixed reciprocal-real space representation 
of the charge density. This allows us to naturally sat- 
isfy the boundary conditions in the different directions. 
Screening functions or other approximations are thus not 
needed. For the direction with the free boundary con- 
ditions, the guidelines described in ref.(0) are followed, 
properly adapted to the surfaces case. 

The charge density in the non-periodic direction is rep- 
resented using interpolating scaling functions, thus avoid- 
ing from the beginning long range interactions between 
super-cells. Like the free BC case, this Poisson solver 
is most efficient when dealing with localized densities 
(in the non-periodic direction). Such densities are for 
instance obtained from electronic structure codes using 
finite differences (fill ) or finite elements (fill ) or also Gaus- 
sians (fl5| ) for the representation of the wave functions. 



in z. Without loss of generality, a function / that lives 
in such a domain can be expanded as 



f(x,y,z) 



e L * L « 'f Px , Py {z) ■ (4) 

Px,P„ 



We indicate with f P:c , Py {z) the one-dimensional function 
associated to the vector p = (p x / L Xl p y /L y ) in the re- 
ciprocal space of the two dimensional surface. Following 
these conventions, the Poisson equation ([TJ become a re- 
lation between the reciprocal space components of V and 
P- 

/+oo 
&z'G{2n\^-z-z')p p ^ Py (z) , (5) 
-oo 

where \p\ 2 = (p x /L x ) 2 + (p y /L y ) 2 , and G(/z; z) is the 
Green function of the one-dimensional Helmholtz equa- 
tion: 



(d 2 -fi 2 ) G(fi;z) = 5(z) 



(6) 



The free BC on the z direction fix the form of the Green 
function: 



{_J_ p -mM ,/ > n 
±\z\ /! = 



(7) 



II. INTERPOLATING SCALING FUNCTIONS 

Interpolating scaling functions (fioh arise in the frame- 
work of wavelet theory They are one-dimensional 
functions, and their three main properties are: 

• The full basis set can be obtained from all the trans- 
lations by a certain grid spacing h of the mother 
function <f> centered at the origin. 



They satisfy the refinement relation 



(x) = J2 h J ^ 2x - 3) 



(3) 



where the hj's are the elements of a filter that char- 
acterizes the wavelet family, and m is the order of 
the scaling function. Eq. ([3]) establishes a relation 
between the scaling functions on a grid with grid 
spacing h and another one with spacing h/2. 

• The mother function <fi is symmetric, with compact 
support from — m to m. It is equal to one at the 
origin and to zero at all other integer points (in 
grid spacing units). The expansion coefficients of 
any function in this basis are just the values of the 
function on the grid. 



III. POISSON EQUATION FOR SURFACES BC 



In numerical calculations continuous charge distribu- 
tions are typically represented by their values on a grid. 
The mixed representation of the charge density given 
above immediately suggests to use a plane wave ex- 
pansion in the periodic directions, which may be eas- 
ily treated with conventional FFT techniques. For the 
non-periodic direction z we will use interpolating scaling 
functions representation. The corresponding continuous 
charge distribution is thus given by 



E E 



x exp 



2 Py— — 

-2wi(^-x 



N z 

E/ PPmPyii* 
j*=0 



i. > (8) 



where h is the grid spacing in the z direction, and 
4>(j) — Sjfi, j £ Z. This mixed representation of the 
charge density in eq.© has an important property. It is 
in a certain sense the most faithful continuous charge dis- 
tribution that can be obtained from values on a grid. The 
main features of the electrostatic potential potential are 
determined by the multipoles and Fourier components of 
the charge distribution. For surface problems the multi- 
poles are the significant quantity along the non-periodic 
direction and the Fourier components are the relevant 
quantities in the two periodic directions. Using the com- 
pleteness relation of plane waves (which holds also on a 
discrete real space grid) and the relation 



Consider a three-dimensional domain, periodic (with 
period L x and L y ) in x and y directions, and non-periodic 



3)x =3 



(9) 
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which is demonstrated in (@), it is easy to show that 



dx dy dz cos ( — ^— ) cos ( — ) z l3 p(x, y, z) 



N- 



4 2^ J 3 P^i^2-,h - 

ei,2 = ±l js=0 

JV, N v N z 

= LxLy h^ E E E >< 

n =0j 2 =0j 3 =0 

xcosl — — | cos I ~ — I J3 3 p{—3i,--z-n, h3z) 



V N x 



>N X JX ' N y 



(10) 



with l\ < N x , £2 < N y and £3 < m, where m is the order 
of the scaling function. This equation shows that the 
discrete multipoles and Fourier components of a charge 
density given on a grid are identical to the multipoles and 
Fourier components of the continuous p given by eq. (|8j). 
For this reason it is not necessary that the input charge 
density for our method is given in the representation of 
eq. ([8]). It can, and it will actually in most cases, be given 
by numerical values on a grid. 

Combining eq. ([5]) with ([5]), the discretized Poisson 
problem thus becomes 



where the quantity (kernel) 



j) = d " G(p; h(j - u))cj>{u) 



(11) 



(12) 



is defined via an integral in the dimensionless variable 
u. Due to the symmetry of <p, the kernel is symmetric 
in the non-periodic direction K(p;j z ) — K(p; —j z ). The 
integral bounds can be restricted from — m to m, thanks 
to the compact support of 4>- 

Once we have calculated the kernel, which will be de- 
scribed below, our numerical procedure is the follow- 
ing. We perform a two-dimensional FFT on our real 
space charge density to obtain the Fourier coefficients 
Ppx,p y \j'z l° r au the periodic planes. Then we have to 
solve eq. (|lll) . Since this equation is a convolution it can 
be calculated by zero-padded FFT's. Finally the poten- 
tial is transformed back from the mixed representation 
to real space to obtain the potential on the grid by an- 
other two-dimensional FFT. Due to the FFT's, the total 
computational cost is 0(N log N). Since all quantities 
are real, the amount of memory and the number of oper- 
ations for the FFT can reduced by using real-to-complcx 
FFT's instead of complex-complex FFT's. 

It remains now to calculate the values of the kernel 
function K(ji;j). The core of the calculation is repre- 
sented by the function 




x \ u -i\(t>(u) 



A > , 
A = 0. 



(13) 



The kernel has the properties K(p;j) = —k(ph;j)/(2p) 
for p > and K(0;j) = K(0;j)/2. A simple numerical 



integration with the trapezoidal rule is inefficient since 
G{p; z) is not smooth in z = while the scaling function 
varies significantly around the integer points. Thanks to 
the compact support of the scaling function, this problem 
can be circumvented with a simple and efficient recursive 
algorithm. We define two functions k^ + > and such 
that K(\;j) = K {+) (\;j) + k(~\X;j), where we have, 
for A > 



X«(A;j) = 



dwe~ A ^'V(u) 



while with A = 



A+) 



(?) 



duu 4>(u) 

JO 

dxu (j>(u) 



1=0,1 



(14) 
(15) 

(16) 
(17) 

(18) 



These objects satisfy recursion relations: 

KW(\;j + 1) = \k^(X;j) ± e^D^(j) 



zj, ±) (j + i) = z( ±) (j)±c e (j), e = o,i, (19) 



where 



D 



Ct(j) 



due ±Xu 4>(u) 



du \ 



(«), £ = 0,1. 



From eq . (fTH - l2Tj) . and the properties 
k(X; 3) = K{\- -j) , KM (A; 0) = K^(X; 0) 



(20) 
(21) 

(22) 



A+) 



(0) 



(0) 



(0) = 4" ) (0) 



the function k(X;j) can be calculated recursively for 
each j e N, by knowing £(+)(A;0) and z[ +) {0), then 
evaluating D^\j) and Q(i) for each value of j. The 
integrals involved can be calculated with high accu- 
racy with a simple higher-order polynomial quadrature. 
They are integral of smooth, well-defined quantities, since 
the interpolating scaling function goes to zero at their 
bounds. Moreover, for values of j lying outside the sup- 
port of <f> we can benefit of a functional relation for cal- 
culating the values of the kernel. The support of a m-th 
order scaling function goes from — m to m, then we have 
Vp > 



K(p; m+p) = e-^ hp K{p- m) 
K(0;m + p) = K(0;m 



p > 



V Zq +) (m) 



(23) 



To summarize, we have found an efficient method for 
evaluating equation (fl"2|) for 3 = 0, • • • , N z and a fixed 
p. Instead of calculating N z + 1 integrals of range 2m, 
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we can obtain the same result by calculating 2 integrals 
of range m and 4m integrals of range 1, with the help of 
relation (|23|) . This will also increase accuracy, since the 
integrands are always smooth functions, which would not 
be the case with a naive approach. 

The accuracy in calculating the integrals can be fur- 
ther improved by using the refinement relation §3§ for 
interpolating scaling functions. For positive A we have 

K(2X;i) = J due- 2X \ u - % \(j){u) 

= ~ J due- x \ u - 2l \(/>(u/2) 

= ^I>/ ^e-^- 2 Mu-j) (24) 

3 

= i£>if(A;2i-j). 

3 

This relation is useful to improve the accuracy in evalu- 
ating the kernel for high A. Since in this case the expo- 
nential function is very sharp, it is better to calculate the 
kernel for lower A and an enlarged domain and then ap- 
ply relation as many times as needed. The relation 
(l23f allows us to enlarge the domain with no additional 
computational cost. With the help of the above described 
properties the computational time for evaluating the ker- 
nel in Fourier space can be considerably optimized, be- 
coming roughly half of the time needed for its application 
on a real space density. 

IV. NUMERICAL RESULTS AND COMPARISON WITH 
OTHER METHODS 

Our new method has an algebraic convergence rate 
of h m , where h is the grid spacing and m is the or- 
der of the interpolating scaling function used. This is 
due to the fact that the charge density along the non- 
periodic direction can be represented with an approxi- 
mation error of the order of h m using interpolating scal- 
ing functions of order m (fl7l ). The error of the plane 
wave representation of the charge density in the peri- 
odic plane has a faster exponential convergence rate and 
is thus not visible in the overall convergence rate. By 
choosing higher order interpolating scaling function we 
can get arbitrarily high algebraic convergence rates. Our 
method was compared with the reciprocal space meth- 
ods by Hockney (@) and Mortensen |j) as implemented 
in the CPMD electronic structure program (16). The 
tests shown in figure [1] are performed with an analytical 
charge distribution that is the Laplacian of V(x, y, z) — 
exp(cos(|^a;) + cos(|f4/)) exp(- -^jj - tan(f-z) 2 ). Its 
behavior along the xy surface is fully periodic, with all 
the reciprocal space components taken into account. The 
exp(— tan(-^-z) 2 ) guarantees a localized behavior in the 
non-periodic direction with the potential going explicitly 
to zero at the borders. This makes also this function 
suitable for comparison with reciprocal space based ap- 
proach. The Gaussian factor is added to suppress high 
frequency components. Tests with other analytical func- 
tions gave comparable accuracies. The reciprocal space 
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FIG. 1 Accuracy comparison between our method with scal- 
ing functions of different orders and the Mortensen solver for 
surface systems as implemented in CPMD. The results of the 
Hockney method are not shown since they are much less pre- 
cise. The h 8 curve is plotted to show the algebraic decrease of 
the precision with respect to the grid space h. The accuracy 
is finally limited by the evaluation of the integral (|13[) . which 
is computed with nearly machine precision. 

Poisson solvers turn out to be much less precise than 
our approach, which explicitly preserves the BC along 
each direction. Moreover, the accuracy shown for the 
Mortensen approach is calculated only for planes that 
lies in the bulk of the non-periodic direction (30% of the 
total volume). Outside of this region, errors blow up, 
which will imply that a very large computational volume 
is needed to obtain accurate results in a sufficiently large 
domain of the non-periodic direction. 

To show that our method genuinely preserves the 
boundary conditions appropriate for surfaces we calcu- 
lated the electrostatic potential for a plane capacitor. For 
this system only the zero-th Fourier components in the 
plane are non- vanishing. Figure[2]shows the results cither 
in the Mortensen/ Hockney reciprocal space methods or 
with our approach. For the plane capacitor, the screening 
function used in the Mortensen approach vanishes, and 
the solution is equal to what we would have obtained with 
a fully periodic boundary conditions. To obtain the good 
"zero electric field" behavior in the borders that we ob- 
tain directly with our method one would have to postpro- 
cess the solution obtained from the Mortensen method, 
by adding to the potential a suitable linear function along 
the non-periodic direction. This is legitimate since a lin- 
ear function is annihilated by the Laplace operator and 
the modified potential is thus also a valid solution of the 
Poisson equation just with different boundary conditions. 
The Hockney method presents a better qualitative be- 
havior, though the results are not accurate. Only with 
our approach we get both accurate and physically sound 
results. 

Table Q] shows the required CPU time for solving the 
Poisson equation on a grid of 128 3 grid points as a func- 
tion of the number of processors on a Cray XT3 parallel 
computer. The parallel version is based on a parallel 
3-dimensional FFT, where the input/output is properly 
distributed/gathered to all the processors. The FFT's 
are performed using a modified version of the algorithm 
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FIG. 2 Electrostatic potential V for a system with two pe- 
riodic planes charged with opposite sign (plane capacitor), 
oriented along the z direction, calculated by different Poisson 
solvers. The values of V are taken in the middle of the x (pe- 
riodic) direction. The position of the positive (black, dashed 
line) and the negative (red, solid) charged plane is schemati- 
cally shown in the figure. The represented solution are, from 
top to bottom, the results from the Mortensen, the Hockney 
and our Poisson solver. 

described in ref.(9) that gives high performances on a 
wide range of computers. 

A package for solving Poisson's equa- 
tion in real space according to the method 
described here can be downloaded from 
http://www.unibas.ch/comphys/comphys/SOFTWARE 

In conclusion, we have presented a method that al- 
lows us to obtain accurate potentials arising from charge 
distributions on surfaces with a O(AlogiV) scaling in a 
mathematically clean way. This method preserves ex- 
plicitly the required boundary conditions, and can easily 
be used for applications inside electronic structure codes 
where the charge density is either given in reciprocal or 
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TABLE I The elapsed time in seconds required on a Cray 
XT3 (based on AMD Opteron processors) to solve Poisson's 
equation with surface BC on a 128 3 grid as a function of 
the number of processors. The time for setting up the Ker- 
nel (around 50% of the total time) is not included. For a 
large number of processors, the communication time needed 
to gather the complete potential to all the processors becomes 
dominant. 



in real space. 
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